Surface Polaron Formation in the Holstein model 
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The effect of a solid-vacuum interface on tlie properties of a strongly coupled electron-phonon system is 
analyzed using dynamical mean-field theory to solve the Holstein model in a semi-infinite cubic lattice. Polaron 
formation is found to occur more easily (i.e., for a weaker electron-phonon coupling) on the surface than in the 
bulk. On the other hand, the metal-insulator transition associated to the binding of polarons takes place at a 
unique critical strength in the bulk and at the surface. 
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I. INTRODUCTION 

Convincing experimental evidence of polaronic behavior 
has been reported in materials such as the high-Tc cuprates 
and manganites. For instance, the transition from the low 
temperature ferromagnetic metallic state to the paramagnetic 
insulating state in manganites is caused by the formation of 
combined structural/magnetic polarons'. Signatures of small 
polarons have been observed in undoped cuprates^'^. 

From a theoretical point of view, polaron formation has 
been intensively studied using a number of approaches. An 
important role to improve our understanding of polaronic phe- 
nomena has been recently played by Dynamical Mean-Field 
Theory (DMFT)'*. DMFT is a powerful non-perturbative tool 
for strongly interacting systems. This technique, which be- 
comes exact in the limit of infinite coordination number, re- 
duces the full lattice many-body problem to a local impurity 
embedded in a self-consistent effective bath of free electrons. 

DMFT studies of the half-filled Holstein model in a Bethe 
lattice with a semi-elliptic free density of states have clarified 
the difference between the polaron crossover, i.e., the con- 
tinuous and progressive entanglement between electrons and 
phonons, and the bipolaronic metal-insulator transition^ '''^. If 
no symmetry breaking is allowed, for small e-ph couplings the 
ground state is metallic with Fermi liquid characteristic. Upon 
increasing the e-ph coupling, the carriers lose mobility, even- 
tually acquiring polaronic character, with a finite lattice distor- 
tion associated to the electron motion. Polaron formation oc- 
curs as a continuous crossover. Once formed, polarons tend to 
attract and form a bound pair in real space, called bipolaron^. 
Within DMFT the bipolaronic binding gives rise to an insu- 
lating state of localized pairs ^, and bipolaron formation gives 
rise to a metal-insulator transition. The pairing transition does 
not coincide with the polaron crossover: Polarons are formed 
before (i.e., for a weaker coupling) the pairing transition oc- 
curs as long as the typical phonon frequency is smaller than 
the electronic energy scales (adiabatic regime). 

On the other hand, fabrication of a variety of hetrostruc- 
tures and interfaces involving cuprates and manganites raises 



the question of whether the electronic behavior at the surface 
or interface is different from the bulk. Several studies have 
been devoted to the case of repulsive electron-electron interac- 
tions (Hubbard model) and to a vacuum-solid interface. Pot- 
thoff and Nolting^, and Liebsch"^ have argued that reduced 
coordination at the surface may enhance correlation effects. 
They also studied the magnetic ordering induced by enhanced 
correlation at the surface. Matzdorf et at. ' ' proposed that fer- 
romagnetic ordering is stabilized at the surface by a lattice 
distortion. Surface ferromagnetism had been also discussed 
in a dynamical mean field theory of Hubbard model by Pot- 
thoff andNolting'-. Helmesefa/. studied the scaling behavior 
of the metallic penetration depth into the Mott insulator near 
the critical Coulomb interaction within the Hubbard model'^. 
Borghi et al. have shown the existence of a dead surface layer 
with exponentially suppressed quasiparticles'"^. 

Here we concentrate on the effect of a vacuum-solid sur- 
face on interacting electron-phonon systems, and we study the 
formation of polarons and the transition to a bipolaronic insu- 
lating state in the semi-infinite Holstein model at half-filling 
and zero temperature on the bipartite simple cubic (sc) lat- 
tice with nearest-neighbor hopping. While the occurrence of 
charge transfer is typical for a system with reduced transla- 
tional symmetry, in our model at half-filling any charge trans- 
fer is excluded by the particle-hole symmetry, leading to a 
homogeneous charge distribution among the layers parallel to 
the surface, and local occupations near the surface do not dif- 
fer from the average filling, {no) — (n) = 1, where a labels 
each layer. 

In addition to the geometrical effect of missing neighbors, 
the surface electronic structure of interacting electron systems 
is also complicated by the fact that the microscopic interac- 
tions in the vicinity of the surface have values which may sig- 
nificantly differ from those in the bulk. A relaxation of the 
surface layer, for example, changes the overlap between the 
one-particle basis states and thus implies a modified hopping 
integral. Within the Holstein model, the parameter modifica- 
tions will be reflected in different values of the surface top- 
most layer hopping integrals and e-ph coupling strengths rel- 
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ative to the bulk ones. In this work we will not consider this 
effect, in order to focus on the more intrinsic effects that can 
give rise to different physics in the surface. 

This paper is organized as follows: In Sec. II the model is 
introduced and the application of DMFT for surface geometry 
is briefly discussed. We mainly characterized the electronic 
and phononic properties by considering the layer quasiparti- 
cle weights, double occupancies and the phonon probabiUty 
distribution function. The corresponding results are discussed 
in Sec. III. Finally Sec. IV is devoted to concluding remarks. 



II. THE MODEL AND METHOD 



The Holstein Hamiltonian is defined by: 



+ gY,{n, - 1) (b\ + b,) + f7o ^ b\h,, (1) 



where Ci^ {^c\a^ ™d bi are, respectively, destruction 
(creation) operators for itinerant electrons with spin a and lo- 
cal vibrons of frequency ilg = 0.2i on site i, Ui is the electron 
density on site i, t stands for the itinerant electrons hopping 
matrix elements between the nearest-neighbor sites, and g de- 
notes the electron-phonon coupling. We fix the energy scale 
by setting t = 1. 

To obtain the ground state properties of this model, we use 
the embedding approach introduced by Ishida and Liebsch'^ 
to extend DMFT to inhomogeneous systems. In this scheme, 
the system is divided into two parts: The surface region which 
includes the first N layers, and the adjacent semi-infinite bulk 
region (substrate) which is coupled to it. Next, we represent 
the effects of the substrate on the surface region by a complex, 
energy-dependent, embedding potential acting on the Hamil- 
tonian matrix of the surface region. The embedding method re- 
quires to consider a relatively small number of surface layers 
and it is therefore a computationally less expensive extension 
of DMFT in the presence of an interface as compared to the 
slab method, in which the inhomogeneous system is simply 
represented as a finite number of layers'^. 

Because of translational symmetry in the plane perpendic- 
ular to the interface, the embedding potential of the substrate 
is diagonal with respect to the two-dimensional wave vector 
k = {kx, ky) and can be expressed as an TV x matrix by 



S(k, iun) = TG(k,iw„)T, 



(2) 



where G(k, icj„) is the Green's function of the substrate de- 
fined by 

G(k,icj„) = [(zc^„ + - e(k) - ^{iu,,)Y\ (3) 

In here, S(iLj„) is the bulk self-energy, which in the frame- 
work of single-site DMFT, is independent of wave vectors, k, 
and LOn are the Matsubara frequencies. We obtain the self- 
energy by performing a standard DMFT calculation for the 



bulk crystal corresponding to the substrate, /i is the chemical 
potential and e(k) is the two-dimensional dispersion relation, 
which includes information about surface geometry. The e(k) 
matrix for the surface cutting a simple cubic lattice along the 
z direction [sc(OOl) surface] takes the following form^: 



e(k) 




te_L(k) 

te||(k) te_L(k) 

fei(k) te||(k) ••• 





(4) 



The intralayer (parallel) hopping and the interlayer (perpen- 
dicular) hopping are specified by te\\ (k) and iex(k), respec- 
tively, and are given by 



eii -2[cos{kx) + cos(fcy)], |e±(k)p = 1. 



(5) 



Finally, T is the hopping matrix between primitive cells of 
substrate and surface region. Since T is non-zero between 
nearest-neighbor layers of substrate and surface region, only 
the surface Green's function'^ of the substrate need to be con- 
sidered in Eq. (2). 

After constructing the embedding potential of the sub- 
strate, S(k, icc)„), by way of a coupled-layer DMFT calcu- 
lation in the surface region the self-energy matrix is deter- 
mined self-consistently. This can be achieved via the fol- 
lowing steps: (i) associating an effective impurity model with 
each layer in the surface region, solving them by using an im- 
purity solver to find the layer-dependent local self-energies, 
Yia{iijJn), and constructing the surface region self-energy ma- 
trix which is diagonal in layer indices (a, (3) with the ele- 
ments, Yiaiiiiuin) = '^a{i^n)5ai3, (ii) Calculating the on-site 
layer-dependent Green's function via the following relation: 



Ga{iuJn) 

E 



, (6) 



{iujn + (U)l - e(k) - S(k, zw„) - S(iw„) / 



where N x N e(k) matrix is given by Eq. (4), (iii) imple- 
menting the DMFT self-consistency relation for each layer, 

Q^{iujn) = [G^^{iijJn) + SQ,(ia'„)] , which determines the 
bath parameters for the new effective impurity model. The 
cycles have to be repeated until self-consistency is achieved. 

We use the exact diagonalization (ED) technique to solve 
the effective impurity model at zero temperature'^, which 
works equally well for any values of the parameters and only 
involves a discretization of the bath hybridization function, 
which is described in terms of a finite and small set of levels 
ris in order to limit the Hilbert space to a workable size. For 
the case of phonon degrees of freedom we considered here, 
the infinite phonon space is also truncated allowing for a max- 
imum number of excited phonons Nph- The typical values we 
considered for the bath levels are ji^ = 8 — 9 and typical max- 
imum number of phonons are Nph = 30 — 50. We tested that 
these numbers, indeed, provide converged results. Moreover, 
the number of surface layers is chosen to be = 5 in all 
calculations. 
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FIG. 1: Quasiparticle weight z of semi-infinite Holstein model for 
simple cubic lattice in the (001) orientation as a function of layer 
index a. Crosses on the vertical axis on the right hand side indicate 
the bulk 2 corresponding to five given- values of g. Lines are drawn 
as a guide to the eye. 



III. RESULTS 

As we anticipated in the introduction, all our calculations 
are performed for the case of uniform parameters. This as- 
sumption, together with the half-filling condition which en- 
forces charge homogeneity, let us to single out the effect of 
the interface and to focus on the purely geometrical aspect of 
the problem. Fig. 1 shows the calculated quasiparticle weight 
Zq, of the semi-infinite Holstein model at T = in the metal- 
lic range as a function of layer index a, where the outermost 
layer corresponds to a = 1. Zq. measures the metallic nature 
of a system, z being one for a non-interacting metal and zero 
for a correlated insulator. In our case z = implies a bipo- 
laronic insulator. The crosses on the vertical axis on the right 
hand side indicate the z values of the bulk metal determined 
by a separate bulk DMFT calculation. For any value of the 
coupling, the quasiparticle weight of the surface layer Zq=i is 
significantly reduced compared to Za=2 and Zq=3 which can 
be understood as the effect of the reduced surface coordina- 
tion number and enhanced effective correlations, in complete 
analogy to the results for repulsive interactions. The evolution 
as a function of the layer index depends instead on the cou- 
pling regime. For weak and moderate coupling, z^ has a non 
monotonic behavior which is damped with increasing distance 
to the surface. For g-values closer to the critical coupling 
strength of the bulk bipolaronic transition, (<7c ~ 0.55) 
the behavior changes qualitatively. Here the layer dependence 
becomes monotonic, and the quasiparticle weight quickly ap- 
proaches its bulk value with increasing a. 

Fig. 2 illustrates the layer-dependent quasiparticle weight 
Za for the first three layers and the bulk quasiparticle weight 
as a function of g. As expected, all the z's monotonically de- 
crease as a function of the e-ph coupling, and they eventually 
vanish. As can be seen in the figure, the differences between 
the and the bulk z diminish with increasing distance from 
the surface and for the third layer, the quasiparticle weight is 
almost indistinguishable from the bulk z on the scale used. 
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FIG. 2: Top panel : Layer-dependent quasiparticle weight for 
the first three layers of the semi-infinite Holstein model with simple 
cubic (001) surface geometry and the bulk quasiparticle weight as a 
function of e-ph coupling strength, g. a. = \ stand for the topmost 
surface layer. The solid line shows z for bulk calculations. The inset 
shows Za (g) in the critical regime. Bottom panel : Layer-dependent 
double occupancy da as a function of g. 



It is crucial to observe that the different Zq, all vanish at the 
same value of g, which also coincides with the bulk critical 
coupUng strength, gc = gcMik- 

If the surface and the bulk were decoupled, the reduced sur- 
face coordination number would tend to drive the surface to 
an insulating phase at a coupling strength lower than the bulk 
critical coupling gc- However, below gc the bulk excitations, 
due to hopping processes between the surface and the bulk can 
induce a quasiparticle peak with a non-zero weight Za=i > 
in the topmost layer and a real surface transition is not found; 
Za=i remain non-zero, although being very small, up to the 
critical coupling for bulk transition, gc- The investigation 
of the imaginary part of the layer-dependent self-energies, 
'Sa(ii^n) at —^0 (not shown) confirms the uniqueness 
of the critical strength gc- In the limit of ^ and for 
all g < gc, the imaginary part of self energy vanishes for all 
layers as it must be for a Fermi liquid. In the coupling con- 
stants close to gc and in the metallic regime, a significant layer 
dependence of Im'Sa{ii^n) for (jJ„ — > with a considerably 
larger slope in the first layer (a = 1) is seen, which reflects the 
enhanced correlation effects at the surface. In the insulating 
state, ImT^a{ioJn) diverges for uun ^ 0. Therefore, there is 
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a unique critical strength gc at which all quasiparticle weight 
functions, Zc((.g), simultaneously approach zero. 

The layer-dependent average double occupancy, = 
{naiUai), is shown in the bottom panel of Fig. 2 as a func- 
tion of g. For small g and all layers, da increases gradually. 
At g = gc it rapidly reaches w 1/2. In the metallic region, 
the double occupancies are increased more rapidly at the top- 
most surface layer as compared with the interior of the sys- 
tem. Again, this is due to the stronger effective e-ph interac- 
tion which results from the narrowing of the non-interactiong 
density of states at the surface. 

We have thus far established that even in the presence of 
a surface, the half-filled Holstein model undergoes a single 
bipolaronic metal-insulator transition, despite the surface is 
less metallic than the bulk for any g < g^- We now discuss 
how the surface influences the local lattice distortions, mea- 
sured by the phonon probability distribution function (PDF), 
P{x) ~ {(f)o\x) {x\4>o) , where |a;)(.T| is the projection operator 
on the subspace for which the phonon displacement x has a 
given value x, and {(po) is the ground state vector. This quan- 
tity can be used to characterize the polaron crossover'^. 

In the absence of e-ph interaction P{x) is a Gaussian cen- 
tered around x = 0. A small e-ph coupling slightly broad- 
ens the distribution which remains centered around a; = 0, 
implying that the coupling is not sufficient to give rise to a 
finite polarization of the lattice. Continuously increasing the 
interaction one eventually obtains a bimodal distribution with 
two identical maxima at x = ±xo. Those maxima are in- 
deed associated with empty and doubly occupied sites, and 
testify the entanglement between the electronic state and the 
lattice distortion, which is precisely the essence of the polaron 
crossover Thus, the appearance of a bimodal shape in P{x) 
is a marker of the polaron crossover^ '^. 

Fig. 3 shows the polaron crossover for our semi-infinite 
Holstein model. For each layer the evolution as a function of 
the coupling follows the pattern we described above: The an- 
harmonicity due to e-ph interaction increases with increasing 
coupling strength leading first to a non-Gaussian and finally 
to a bimodal PDF at all g > gpoi . This behavior signals the 
appearance of static distortions, even if we are neglecting any 
ordering between them. The strongest differences with respect 
to the bulk PDF are found for the top layer (a = 1) PDF. The 
layer PDFs converge to the bulk PDF with increasing distance 
to the surface. Beyond the third layer the PDF is essentially 
identical to its bulk behavior. It is apparent from the data of 
Fig. 3 that the PDF at the topmost (surface) layer becomes 
bimodal at lower values of the coupling strength with respect 
to the internal layers, the surface can display polaronic dis- 
tortions while the bulk is still undistorted (even if the local 
vibrations are strongly anharmonic). 



IV. CONCLUDING REMARKS 

We have investigated polaron formation and transition to 
the bipolaronic insulating state at solid- vacuum surface at zero 
temperature in the framework of the semi-infinite Holstein 
model at half-filling. Using the embedding approach to extend 
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FIG. 3: Phonon probability distribution function for the first and 
second surface layers of a semi-infinite sc(OOl) Holstein model at 
half-filhng. The various curves refer to different values of electron- 
phonon coupling strength, g. Upon increasing the e-ph coupling, 
a smooth crossover occurs between a unimodal distribution and a 
bimodal distribution. However, the polaron crossover (onset of bi- 
modality) occurs at different values of g for the first and second lay- 
ers. Top panel shows that at the topmost surface layer, polaron for- 
mation takes place at g ~ 0.51, whereas at the second layer, it takes 
place at g ~ 0.53 (bottom panel). All other layers behave just like 
the second layer. 



dynamical mean-field theory to layered systems, it is found 
that the bipolaronic insulating state occurs simultaneously at 
the surface and in the bulk, and it takes place exactly at the 
same critical coupling strength gc,buik as for the infinitely ex- 
tended system, gc,buik = gc- When the system is metallic 
the topmost layer quasiparticle weight zi is smaller than the 
bulk value z^uik, since a reduced surface coordination num- 
ber implies a stronger effective correlation effects. Fixing 
the coupling at values quite smaller than g^,, the quasiparti- 
cle weight is an oscillating function of the layer index. As 
the distance from the surface increases, these oscillations fade 
away. For couplings close to the metal-insulator transition Za 
instead monotonically increases by approaching the bulk. On 
the other hand, the polaron crossover occurs more easily at 
the surface with respect to the bulk. There is therefore a fi- 
nite window of e-ph coupling in which the surface presents 
polaronic distortions, while the bulk has no distortions. As 
we already mentioned, this difference is not able to support a 
metallic bulk coexisting with an insulating surface. 



5 



Surface effects are expected to be the more pronounced 
the larger is the number of missing neighbors in the topmost 
layer. As we move from the sc(OOl) to the sc(01 1) and to the 
sc(l 1 1) surface geometries, the surface coordination numbers 
decrease from nc"°^'' = 5 to ni"^^-* = 4 to ni^^'^^ = 3, respec- 
tively. Therefore, we expect to observe a narrowed topmost 
layer free density of state which results in an enhanced ra- 
tio between the e-ph coupling strength and the effective band 
width. Consequently, the e-ph interaction tends to be stronger 
at the surface. Clearly, according to this argument we expect 
the difference between the two coupling strengths for the po- 
laron formation at the surface and in the bulk would become 



larger 

The present study has been restricted to uniform model pa- 
rameters. This leaves several open questions like the possi- 
bility of coexisting different surface and bulk phases, if the 
model parameters at the vicinity of the surface are modified 
(See the comment made in this respect in Sec. I). 
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